Source code and data found at https://github.com/maRce10/budgie_inbre_stress_vocal_output

Load packages

## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2",
    "kableExtra", "knitr", "formatR", "MASS", "sp", "GGally", "brms",
    "lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes",
    "cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")

source("~/Dropbox/R_package_testing/sketchy/R/load_packages.R")
load_packages(x)

source("~/Dropbox/R_package_testing/brmsish/R/html_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/read_summary.R")

Functions and global parameters

opts_knit$set(root.dir = "..")

# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE,
    tidy = TRUE)

read_excel_df <- function(...) data.frame(read_excel(...))


# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")


standard_error <- function(x) sd(x)/sqrt(length(x))

cols <- viridis(10, alpha = 0.7)

col_pointrange <- cols[7]
# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")

# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd",
    "freq.median", "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent",
    "time.ent", "entropy", "meandom", "mindom", "maxdom", "dfrange",
    "modindx", "meanpeakf")]

sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])

sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year,
    sep = "-")

sels$date <- as.Date(sels$date, format = "%d-%b-%Y")

# acoustic measurements
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")

indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")

# group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID),
# function(x) { X <- group_ovlp[group_ovlp$ID == x, ]
# X$overlap.to.group <- X$overlap.to.group -
# X$overlap.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] X$distance.to.group <-
# X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))

group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

Counts per individual

df <- as.data.frame(table(sels$ID))

names(df) <- c("ID", "Sample_size")

df <- df[order(df$Sample_size, decreasing = FALSE), ]

kb <- kable(df, row.names = FALSE)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover",
    "condensed", "responsive"))

print(kb)
ID Sample_size
132YGMM 6
125YGMM 12
160YGHM 16
310YGHM 21
393YGHM 25
303YGHM 37
398WBHM 41
365YGHM 44
399YGLM 46
300YGHM 47
270YGHM 64
407YGHM 97
386WBHM 100
377WWLM 108
367WWNM 119
371YYLM 149
397YGHM 175
378YYLM 195
404WBHM 196
258YGHM 223
389WWLM 230
262YGHM 306
400YGHM 360
388YGLM 377
327YYLM 404
396YBHM 455
403WBHM 566
356WBHM 761
361YGHM 767
402YGHM 849
395WBHM 854
360YGHM 900
390WBHM 935
405YBHM 975
385YBHM 1256
394YBHM 1693

Add metadata

metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")

# head(metadat)

sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"

# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID,
# sels$ID)

sels$treatment <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]

})


sels$treatment.room <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]

})


sels$round <- sapply(1:nrow(sels), function(x) {

    metadat$Round[metadat$Bird.ID == sels$ID[x]][1]

})

sels$source.room <- sapply(1:nrow(sels), function(x) {

    metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]

})

sels$record.group <- sapply(1:nrow(sels), function(x) {

    metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]

})


# add week
out <- lapply(unique(sels$round), function(x) {

    Y <- sels[sels$round == x, ]

    min_date <- min(Y$date)
    week_limits <- min_date + seq(0, 100, by = 7)

    Y$week <- NA
    for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i -
        1] & Y$date < week_limits[i]] <- i - 1

    return(Y)
})

sels <- do.call(rbind, out)


sels$cort.baseline <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$cort.stress <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})

sels$stress.response <- sels$cort.stress  #- sels$cort.baseline

sels$weight <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$breath.count <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


# remove week 5
sels <- sels[sels$week != 5, ]
agg_dat <- aggregate(selec ~ ID + week, data = sels, length)

# compare to week 1 agg_dat$call.count <-
# sapply(1:nrow(agg_dat), function(x) { baseline <-
# agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]]
# if (length(baseline) > 0) change <- agg_dat$selec[x] -
# baseline else change <- agg_dat$selec[x] return(change) } )

# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])

agg_dat$selec <- NULL

agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID ==
        agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$stress.response[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$stress.response[sels$week == agg_dat$week[x] &
        sels$ID == agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})


agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID ==
        agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID ==
        agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {

    baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
    current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID ==
        agg_dat$ID[x]]

    if (length(baseline) > 0 & length(current) > 0)
        change <- mean(current) - mean(baseline) else change <- NA

    return(change)
})

agg_dat$acoustic.diversity <- sapply(1:nrow(agg_dat), function(x) {

    area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] &
        areas_by_week$week == agg_dat$week[x]]

    if (length(area) < 1)
        area <- NA

    return(area)
})

agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {

    distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID ==
        agg_dat$ID[x] & indiv_ovlp$week == agg_dat$week[x]]

    if (length(distance) < 1)
        distance <- NA

    return(distance)
})

agg_dat$acustic.plasticity <- sapply(1:nrow(agg_dat), function(x) {

    overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
        indiv_ovlp$week == agg_dat$week[x]]

    plasticity <- 1 - overlap

    if (length(plasticity) < 1)
        plasticity <- NA

    return(plasticity)
})

agg_dat$acoustic.convergence <- sapply(1:nrow(agg_dat), function(x) {

    overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] &
        group_ovlp$week == agg_dat$week[x]]

    if (length(overlap) < 1)
        overlap <- NA

    return(overlap)
})


agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
    agg_dat$ID[x]]))

agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
    agg_dat$ID[x]]))

Physiological parameters

Barplot and effect sizes graph

physio_models <- readRDS("./data/processed/physiological_response_models.RDS")

breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count",
    "D14.Breath.Count", "D21.Breath.Count", "D28.Breath.Count")])

weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.",
    "D14.Bird.Weight..g.", "D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])

cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress",
    "D14.CORT.Stress", "D21.CORT.Stress", "D28.CORT.Stress")])

cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline",
    "D14.CORT.Baseline", "D21.CORT.Baseline", "D28.CORT.Baseline")])


stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round",
    "Treatment.Room")], week = breath.count$ind, breath.count = breath.count$values,
    weight = weight$values, cort.stress = cort.stress$values, cort.baseline = cort.baseline$values,
    stress.response = cort.stress$values - cort.baseline$values)

# head(stress)

stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."),
    "[[", 1), levels = c("D3", "D7", "D14", "D21", "D28"))

stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)

stress$treatment <- factor(stress$Treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

# remove 5th week
stress <- stress[stress$week != "D28", ]


stress_l <- lapply(stress$Bird.ID, function(x) {
    X <- stress[stress$Bird.ID == x, ]

    X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
    X$weight <- X$weight - X$weight[X$week == "D3"]
    X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
    X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week ==
        "D3"]
    X$stress.response <- X$stress.response  #- X$stress.response[X$week == 'D3']

    return(X)
})

stress <- do.call(rbind, stress_l)

agg_stress <- aggregate(cbind(breath.count, weight, stress.response,
    cort.baseline) ~ week + treatment, stress, mean)
agg_stress_se <- aggregate(cbind(breath.count, weight, stress.response,
    cort.baseline) ~ week + treatment, stress, standard_error)

names(agg_stress_se) <- paste(names(agg_stress_se), ".se", sep = "")

agg_stress <- cbind(agg_stress, agg_stress_se[, 3:6])

agg_stress$Week <- 1:4

bs <- 10

gg_breath.count <- ggplot(data = agg_stress, aes(x = Week, y = breath.count,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = breath.count -
    breath.count.se, ymax = breath.count + breath.count.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Mean change in breath\nrate (breaths/min)",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

gg_weight <- ggplot(data = agg_stress, aes(x = Week, y = weight, fill = treatment)) +
    geom_bar(stat = "identity") + geom_errorbar(aes(ymin = weight -
    weight.se, ymax = weight + weight.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
    labs(y = "Mean change in \nweight (grams)", x = "Week") + theme_classic(base_size = bs) +
    theme(legend.position = "none")


gg_cort.baseline <- ggplot(data = agg_stress, aes(x = Week, y = cort.baseline,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = cort.baseline -
    cort.baseline.se, ymax = cort.baseline + cort.baseline.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Mean change in\nbaseline CORT (ng/mL)",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")


gg_stress.response <- ggplot(data = agg_stress, aes(x = Week, y = stress.response,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = stress.response -
    stress.response.se, ymax = stress.response + stress.response.se),
    width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in stress response \nCORT (ng/mL)",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")


gg_coeffs_physio <- lapply(physio_models, function(x) {

    vars <- grep("b_", posterior::variables(x), value = TRUE)
    draws <- posterior::as_draws_array(x, variable = vars)

    coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
        0.975), robust = TRUE)

    coef_table$predictor <- rownames(coef_table)
    coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
    coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
    coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
    coef_table <- coef_table[coef_table$predictor != "Intercept",
        ]

    gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
        y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
        col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
        xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
        ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
            plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
            strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
        ggplot2::labs(x = "Effect size", y = "")

    return(gg_coef)
})

cowplot::plot_grid(gg_breath.count, gg_weight, gg_cort.baseline, gg_stress.response,
    gg_coeffs_physio[[grep("breath", names(gg_coeffs_physio))]] +
        theme_classic(base_size = bs), gg_coeffs_physio[[grep("weight",
        names(gg_coeffs_physio))]] + theme_classic(base_size = bs),
    gg_coeffs_physio[[grep("baseline", names(gg_coeffs_physio))]] +
        theme_classic(base_size = bs), gg_coeffs_physio[[grep("response",
        names(gg_coeffs_physio))]] + theme_classic(base_size = bs),
    nrow = 2, rel_heights = c(1.8, 1))

# try bs = 20 for saving plots

# cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_physiology_70dpi.jpeg',
# width = 25, height = 9) cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_physiology_300dpi.jpeg',
# dpi = 300, width = 25, height = 9)

Stats

Models: Predicted physio measure ~ treatment + week (continuous) + IndRandom

Variables (Difference from Week 1): weight, BR, baseline CORT, Stress CORT, Stress Response

responses <- c("baseline.CORT", "stress.response", "stress.CORT",
    "weight", "breath.rate")

predictors <- c("~ treatment + week + (1|ID) + (1|round)")

formulas <- expand.grid(responses = responses, predictors = predictors,
    stringsAsFactors = FALSE)

vars_to_scale <- c(responses, "week")

# remove week 1
sub_agg_dat <- agg_dat[agg_dat$week != 1, ]

for (i in vars_to_scale) sub_agg_dat[, vars_to_scale] <- scale(sub_agg_dat[,
    vars_to_scale])

physio_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) ==
        formulas$responses[x]]), ]
    sub_dat

    mod <- brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 20000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 4, prior = c(prior(normal(0, 5), "b"), prior(normal(0,
            10), "Intercept"), prior(student_t(3, 0, 10), "sd"), prior(student_t(3,
            0, 10), "sigma")))

    return(mod)
})


names(physio_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(physio_models, "./data/processed/physiological_response_models.RDS")
physio_models <- readRDS("./data/processed/physiological_response_models.RDS")

ggl <- list()
for (x in 1:length(physio_models)) {

    ggl[[length(ggl) + 1]] <- html_summary(physio_models[[x]], gsub.pattern = "b_treatment|b_",
        gsub.replacement = "", highlight = FALSE, remove.intercepts = TRUE,
        model.name = names(physio_models)[x])
    print(ggl[[length(ggl)]])
}

baseline.CORT ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) baseline.CORT ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 406 0 12756.08 17437.15 614553873
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress 0.760 0.108 1.420 1 13169.12 17881.30
MediumStress 0.180 -0.506 0.868 1 12756.08 17437.15
week -0.001 -0.142 0.141 1 27598.67 26244.08

stress.response ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) stress.response ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 2129 0 1831.829 3145.311 1852853265
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress -0.204 -0.889 0.465 1.003 1831.829 3145.311
MediumStress 0.097 -0.667 0.839 1.002 2582.692 10967.477
week -0.152 -0.284 -0.020 1.003 6379.514 15589.581

stress.CORT ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) stress.CORT ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 1015 0 7696.824 8195.241 672971911
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress -0.219 -0.905 0.452 1.001 8677.295 16469.149
MediumStress 0.084 -0.684 0.832 1 7696.824 8195.241
week -0.154 -0.284 -0.022 1 25064.631 25062.987

weight ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) weight ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 1142 0 6459.84 11029.72 351999546
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress -0.319 -0.971 0.328 1.001 10517.09 16418.60
MediumStress -0.120 -0.850 0.598 1.002 6459.84 11029.72
week -0.345 -0.478 -0.211 1 14142.17 16891.28

breath.rate ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) breath.rate ~ treatment + week + (1 | ID) + (1 | round) 20000 4 1 10000 732 0 17760.6 14209.36 1444112821
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress 0.147 -0.156 0.449 1 17760.60 14209.36
MediumStress 0.060 -0.264 0.389 1 18526.45 24064.77
week -0.784 -0.898 -0.670 1 24048.27 26983.78

names(ggl) <- names(physio_models)

 

Takeaways

  • Breath rate decreases gradually with time across after the first week

  • Stress response is higher in “high stress” birds compared to first week

 

Acoustic space projection

t-SNE

scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
    "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent",
    "entropy", "meandom", "mindom", "maxdom", "dfrange", "modindx",
    "meanpeakf")])

tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE,
    max_iter = 5000)

saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")

Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")

sels <- data.frame(sels, Y)

sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) +
    geom_point() + labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) +
    theme_classic(base_size = 25) + guides(colour = guide_legend(override.aes = list(size = 10)))

Behavioral parameters

Barplot and effect sizes pgrah

behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")

agg_call.count <- aggregate(cbind(call.count, acoustic.convergence) ~
    week + treatment, agg_dat, mean)

agg_behav <- aggregate(cbind(acoustic.diversity, acustic.plasticity) ~
    week + treatment, agg_dat, mean)

agg_call.count_se <- aggregate(cbind(call.count, acoustic.convergence) ~
    week + treatment, agg_dat, standard_error)

agg_behav_se <- aggregate(cbind(acoustic.diversity, acustic.plasticity) ~
    week + treatment, agg_dat, standard_error)

agg_behav_se <- merge(agg_call.count_se, agg_behav_se, all = TRUE)

names(agg_behav_se) <- paste(names(agg_behav_se), ".se", sep = "")

agg_behav <- merge(agg_call.count, agg_behav, all = TRUE)

agg_behav <- cbind(agg_behav, agg_behav_se[, 3:6])

bs <- 10

agg_behav$treatment <- factor(agg_behav$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

gg_call.count <- ggplot(data = agg_behav, aes(x = week, y = call.count,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = call.count -
    call.count.se, ymax = call.count + call.count.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Vocal output", x = "Week") +
    theme_classic(base_size = bs) + theme(legend.position = "none")

gg_acoustic.diversity <- ggplot(data = agg_behav, aes(x = week, y = acoustic.diversity,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.diversity -
    acoustic.diversity.se, ymax = acoustic.diversity + acoustic.diversity.se),
    width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Change in vocal diversity",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

gg_acustic.plasticity <- ggplot(data = agg_behav, aes(x = week, y = acustic.plasticity,
    fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acustic.plasticity -
    acustic.plasticity.se, ymax = acustic.plasticity + acustic.plasticity.se),
    width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal plasticity",
    x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")

gg_acoustic.convergence <- ggplot(data = agg_behav, aes(x = week,
    y = acoustic.convergence, fill = treatment)) + geom_bar(stat = "identity") +
    geom_errorbar(aes(ymin = acoustic.convergence - acoustic.convergence.se,
        ymax = acoustic.convergence + acoustic.convergence.se), width = 0.1) +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
    ncol = 3, scale = "fixed") + labs(y = "Vocal convergence", x = "Week") +
    theme_classic(base_size = bs) + theme(legend.position = "none")

gg_coeffs_behav <- lapply(behav_models, function(x) {

    vars <- grep("b_", posterior::variables(x), value = TRUE)
    draws <- posterior::as_draws_array(x, variable = vars)

    coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
        0.975), robust = TRUE)

    coef_table$predictor <- rownames(coef_table)
    coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
    coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
    coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
    coef_table <- coef_table[coef_table$predictor != "Intercept",
        ]

    gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
        y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
        col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
        xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
        ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
            plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
            strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
        ggplot2::labs(x = "Effect size", y = "")

    return(gg_coef)
})

cowplot::plot_grid(gg_call.count, gg_acoustic.diversity, gg_acustic.plasticity,
    gg_acoustic.convergence, gg_coeffs_behav[[grep("count", names(gg_coeffs_behav))]] +
        theme_classic(base_size = bs), gg_coeffs_behav[[grep("diversity",
        names(gg_coeffs_behav))]] + theme_classic(base_size = bs),
    gg_coeffs_behav[[grep("plasticity", names(gg_coeffs_behav))]] +
        theme_classic(base_size = bs), gg_coeffs_behav[[grep("convergence",
        names(gg_coeffs_behav))]] + theme_classic(base_size = bs),
    nrow = 2, rel_heights = c(1.8, 1))

# try bs = 20 for saving plots

# cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_behavior_70dpi.jpeg', width
# = 25, height = 9) # cowplot::ggsave2(filename =
# './output/bar_graphs_and_estimates_behavior_300dpi.jpeg', dpi
# = 300, width = 25, height = 9)

Stats

Model: Predicted behavior ~ treatment + week (continuous) + IndRandom

Variables: # calls, Distance moved from self in first week, Overlap to original acoustic space, Match to group repertoire, Maybe overall size of acoustic space

Do as comparison to week one using rarefacted calls and kernel density

responses <- c("call.count", "acoustic.diversity", "acustic.plasticity",
    "acoustic.convergence")

predictors <- c("~ treatment + week + (1|ID) + (1|round)")

formulas <- expand.grid(responses = responses, predictors = predictors,
    stringsAsFactors = FALSE)

vars_to_scale <- c(responses, "week")


for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[,
    vars_to_scale])

behav_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
        ]

    # remove week 1
    if (!grepl("count|group", formulas$responses[x]))
        sub_dat <- sub_dat[sub_dat$week != 1, ]

    mod <- brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 50000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9,
            max_treedepth = 15), chains = 4, prior = c(prior(normal(0,
            5), "b"), prior(normal(0, 10), "Intercept"), prior(student_t(3,
            0, 10), "sd"), prior(student_t(3, 0, 10), "sigma")))

    return(mod)
})


names(behav_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(behav_models, "./data/processed/behavioral_response_models.RDS")
behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")

ggl <- list()
for (x in 1:length(behav_models)) {

    ggl[[length(ggl) + 1]] <- html_summary(behav_models[[x]], gsub.pattern = "b_treatment|b_",
        gsub.replacement = "", highlight = FALSE, remove.intercepts = TRUE,
        model.name = names(behav_models)[x])
    print(ggl[[length(ggl)]])
}

call.count ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) call.count ~ treatment + week + (1 | ID) + (1 | round) 50000 4 1 25000 1455 0 21916.29 19215.21 2086259884
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress 0.784 0.205 1.357 1 21916.29 19215.21
MediumStress 0.394 -0.183 0.971 1 29122.91 29122.87
week -0.059 -0.203 0.086 1 48490.61 55075.00

acoustic.diversity ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) acoustic.diversity ~ treatment + week + (1 | ID) + (1 | round) 50000 4 1 25000 2953 0 6138.296 3258.34 1000749341
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress -0.208 -1.070 0.657 1.001 6138.296 3258.34
MediumStress -0.699 -1.559 0.158 1.001 8701.168 12053.83
week 0.036 -0.081 0.152 1 27120.995 32683.11

acustic.plasticity ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) acustic.plasticity ~ treatment + week + (1 | ID) + (1 | round) 50000 4 1 25000 7397 0 1025.75 483.942 803966074
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress -0.904 -1.731 -0.062 1.003 4966.382 30252.390
MediumStress -0.232 -1.130 0.678 1.005 1025.750 483.942
week 0.171 0.034 0.306 1.002 2427.679 26782.955

acoustic.convergence ~ treatment + week + (1|ID) + (1|round)

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) acoustic.convergence ~ treatment + week + (1 | ID) + (1 | round) 50000 4 1 25000 8224 0 150.4 57.174 1972081133
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
HighStress 0.335 -0.307 1.003 1.006 1086.629 21463.246
MediumStress 0.372 -0.288 1.045 1.003 12829.194 16680.338
week -0.228 -0.406 -0.058 1.017 150.400 57.174

names(ggl) <- names(behav_models)

 

Takeaways

  • Higher vocal output in “high stress” birds compared to control

  • Higher acoustic overlap to themselves in week 1 for “high stress” birds compared to control

  • Decrease in overlap to themselves in week 1 across time

 


Session information

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] PhenotypeSpace_0.1.0 warbleR_1.1.26       NatureSounds_1.0.4  
##  [4] seewave_2.1.8        tuneR_1.3.3.1        ggridges_0.5.3      
##  [7] posterior_1.2.0      ggrepel_0.9.1        cowplot_1.1.1       
## [10] tidybayes_3.0.2      modelr_0.1.8         tidyr_1.2.0         
## [13] forcats_0.5.1        purrr_0.3.4          dplyr_1.0.8         
## [16] lme4_1.1-28          Matrix_1.3-4         brms_2.16.3         
## [19] Rcpp_1.0.8           GGally_2.1.2         sp_1.4-6            
## [22] MASS_7.3-54          formatR_1.11         knitr_1.37          
## [25] kableExtra_1.3.4     ggplot2_3.3.5        viridis_0.6.2       
## [28] viridisLite_0.4.0    pbapply_1.5-0        readxl_1.3.1        
## [31] lubridate_1.8.0      remotes_2.4.2       
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2            tidyselect_1.1.2      fftw_1.0-6.1         
##   [4] htmlwidgets_1.5.4     grid_4.1.0            munsell_0.5.0        
##   [7] codetools_0.2-18      DT_0.20               miniUI_0.1.1.1       
##  [10] withr_2.4.3           spatstat.random_2.1-0 Brobdingnag_1.2-7    
##  [13] colorspace_2.0-3      highr_0.9             uuid_1.1-0           
##  [16] rstudioapi_0.13       stats4_4.1.0          dtw_1.22-3           
##  [19] tensor_1.5            bayesplot_1.8.1       labeling_0.4.2       
##  [22] emmeans_1.8.1-1       rstan_2.21.3          polyclip_1.10-0      
##  [25] farver_2.1.0          bridgesampling_1.1-2  rprojroot_2.0.2      
##  [28] coda_0.19-4           vctrs_0.3.8           generics_0.1.2       
##  [31] TH.data_1.1-0         xfun_0.29             R6_2.5.1             
##  [34] markdown_1.1          gamm4_0.2-6           projpred_2.0.2       
##  [37] bitops_1.0-7          spatstat.utils_2.3-0  reshape_0.8.8        
##  [40] assertthat_0.2.1      promises_1.2.0.1      scales_1.1.1         
##  [43] multcomp_1.4-17       gtable_0.3.0          goftest_1.2-3        
##  [46] processx_3.5.2        xaringanExtra_0.7.0   sandwich_3.0-1       
##  [49] rlang_1.0.1           systemfonts_1.0.4     splines_4.1.0        
##  [52] spatstat.geom_2.3-2   broom_0.7.12          checkmate_2.0.0      
##  [55] inline_0.3.19         yaml_2.3.5            reshape2_1.4.4       
##  [58] abind_1.4-5           threejs_0.3.3         crosstalk_1.2.0      
##  [61] backports_1.4.1       httpuv_1.6.5          rsconnect_0.8.25     
##  [64] tensorA_0.36.2        tools_4.1.0           ellipsis_0.3.2       
##  [67] spatstat.core_2.4-0   raster_3.5-15         jquerylib_0.1.4      
##  [70] RColorBrewer_1.1-2    proxy_0.4-26          plyr_1.8.6           
##  [73] base64enc_0.1-3       RCurl_1.98-1.6        ps_1.6.0             
##  [76] prettyunits_1.1.1     rpart_4.1-15          deldir_1.0-6         
##  [79] zoo_1.8-9             cluster_2.1.2         magrittr_2.0.2       
##  [82] ggdist_3.1.0          colourpicker_1.1.1    mvtnorm_1.1-3        
##  [85] matrixStats_0.61.0    shinyjs_2.1.0         mime_0.12            
##  [88] evaluate_0.15         arrayhelpers_1.1-0    xtable_1.8-4         
##  [91] shinystan_2.5.0       gridExtra_2.3         rstantools_2.1.1     
##  [94] compiler_4.1.0        tibble_3.1.6          crayon_1.5.0         
##  [97] minqa_1.2.4           StanHeaders_2.21.0-7  htmltools_0.5.2      
## [100] mgcv_1.8-36           later_1.3.0           RcppParallel_5.1.5   
## [103] DBI_1.1.1             boot_1.3-28           permute_0.9-7        
## [106] cli_3.2.0             parallel_4.1.0        igraph_1.2.11        
## [109] pkgconfig_2.0.3       signal_0.7-7          terra_1.5-21         
## [112] spatstat.sparse_2.1-0 xml2_1.3.3            svUnit_1.0.6         
## [115] dygraphs_1.1.1.6      svglite_2.1.0         bslib_0.3.1          
## [118] webshot_0.5.2         estimability_1.4.1    rvest_1.0.2          
## [121] stringr_1.4.0         distributional_0.3.0  callr_3.7.0          
## [124] digest_0.6.29         vegan_2.5-7           spatstat.data_2.1-2  
## [127] rmarkdown_2.11        cellranger_1.1.0      shiny_1.7.1          
## [130] gtools_3.9.2          rjson_0.2.21          nloptr_2.0.0         
## [133] lifecycle_1.0.1       nlme_3.1-152          jsonlite_1.8.0       
## [136] fansi_1.0.2           pillar_1.7.0          lattice_0.20-44      
## [139] loo_2.4.1             fastmap_1.1.0         httr_1.4.2           
## [142] pkgbuild_1.3.1        survival_3.2-11       glue_1.6.1           
## [145] xts_0.12.1            shinythemes_1.2.0     stringi_1.7.6        
## [148] sass_0.4.0